module test
    implicit none
contains

    pure elemental subroutine swap(a,b)
        real(kind=8),intent(inout)::a
        real(kind=8),intent(inout)::b
        real(kind=8)::tmp
        tmp=a;a=b;b=tmp
    end subroutine swap

    subroutine solve(a,b,n,m)
        ! n*n * n*m = n*n
        real(kind=8),parameter::eps=1.d-15
        integer,intent(in)::n
        integer,intent(in)::m
        real(kind=8),intent(inout)::a(n,n)
        real(kind=8),intent(inout)::b(n,m)
        real(kind=8)::temp,temp2
        integer::i,j,k,major
        ! column major method
        do i=1,n-1
            ! 寻找第i列主元 find major position
            major=maxloc(abs(a(i:n,i)),dim=1)+i-1
            !数据交换,cache miss
            if(abs(a(major,i))<eps)error stop "Singular Matrix"
            if(major/=i)then
                call swap(a(i,i:n),a(major,i:n))
                !a([i,major],i:n)=a([major,i],i:n)
                !a(i:major:(major-i),i:n)=a(major:i:i-major,i:n)
                !b([i,major],1:n)=b([major,i],1:n)
                !b(i:major:(major-i),1:n)=b(major:i:i-major,1:n)
                call swap(b(i,:),b(major,:))
            end if
            temp=1.d0/a(i,i)
            !计算
            do j=i+1,n
                a(i,j)=a(i,j)*temp!归一化
                a(i+1:n,j)=a(i+1:n,j)-a(i+1:n,i)*a(i,j)
            end do
            do j=1,m
                b(i,j)=b(i,j)*temp!归一化
                b(i+1:n,j)=b(i+1:n,j)-a(i+1:n,i)*b(i,j)
            end do
            !第i列
            a(i,i)=1.d0
            a(i+1:n,i)=0.d0
        end do
        temp=1.d0/a(n,n)
        a(n,n)=1.d0
        ! 回代
        do j=1,m
            !归一化
            b(n,j)=b(n,j)*temp
            do k=1,n-1
                temp2=b(n-k+1,j)
                b(1:n-k,j)=b(1:n-k,j)-temp2*a(1:n-k,n-k+1)
            end do
        end do
    end subroutine solve
    subroutine gauss (a,n)       ! invert matrix by gauss method
        ! --------------------------------------------------------------------
        implicit none
        integer,intent(in):: n
        real(8),intent(inout):: a(n,n)
        ! - - - local variables - - -
        real(8) :: b(n,n), c, d, temp(n)
        integer :: i, j, k, m, ipiv(n)
        ! - - - - - - - - - - - - - -
        b = a
        ipiv = [(i, i = 1, n)]
        do k = 1,n
            m = maxloc(abs(b(k:n,k)),dim=1)+k-1
            if(abs(b(m,k))<1.d-14) then
                write(*,*) 'Singular Matrix'
                return
            end if
            if (m /= k) then
                ipiv([m,k]) = ipiv([k,m])
                b([m,k],:) = b([k,m],:)
            end if
            d = 1/b(k,k)
            temp = b(:,k)
            do j = 1, n
                c = b(k,j)*d
                b(:,j) = b(:,j)-temp*c
                b(k,j) = c
            end do
            b(:,k) = temp*(-d)
            b(k,k) = d
        end do
        a(:,ipiv) = b
    end subroutine gauss
end module test

program main
    use test
    implicit none
    integer,parameter::n=700,m=700
    real(kind=8)::a(n,n),a0(n,n)
    real(kind=8)::b(n,m),x(n,m)
    integer::ipiv(n),info
    integer::time1,time2
    integer::time_rate
    call random_number(a)
    a=matmul(transpose(a),a)
    call random_number(b)
    !a=reshape([1,2,3,3,2,1,4,4,1],shape=shape(a))
    !b=reshape([1,2,3],shape=shape(b))

    a0=a
    x=b
    call system_clock(time1,time_rate)
    call solve(a, x, n, m)
    call system_clock(time2,time_rate)
    write(*,*)(time2-time1)/real(time_rate,kind=8)
    write(*,*)sum(abs(matmul(a0,x)-b))/(n*m)

    a=a0
    x=b
    call system_clock(time1,time_rate)
    call dgesv(n,m,a,n,ipiv,x,n,info)
    call system_clock(time2,time_rate)
    write(*,*)(time2-time1)/real(time_rate,kind=8)
    write(*,*)sum(abs(matmul(a0,x)-b))/(n*m)

end program main

